
## Helper functions for the analysis steps of
##
## Christopher Adolph, Kenya Amano, Bree Bang-Jensen, 
## Nancy Fullman, Beatrice Magistro, Grace Reinke, Rachel Castellano,## Megan Erickson, and John Wilkerson, "The Pandemic Policy U-Turn:## The role of partisanship, public health, and race in decisions
## to ease COVID-19 social distancing policies in the U.S.", 
## forthcoming in Perspectives on Politics.
##
## Christopher Adolph
## 18 July 2021     faculty.washington.edu/cadolph



## Compute hazard ratios & CIs for continuous covariates from fitted CPH model
simCoxPH <- function(sims, before, after, baseline=0, ci=0.95) {
    crossRR <- sims^(after - before)
    res <- c(mean(crossRR), quantile(crossRR, probs=c((1-ci)/2, (1+ci)/2)))
    res
}



## Plot hazard ratios from CPH models using tile, including special 
## options for robustness checks to show range of results
plotCovidCPH <- function(file, res1, res2=NULL, scenname, xNum, xpre1, xpost1,
                         xpre2=NULL, xpost2=NULL,
                         sublabel1, sublabel2=NULL,
                         limits=c(0.6, 5),
                         at=c(0.75, 1, 1.5, 2, 4),
                         labsX=c("0.75x", "1x", "1.5x", "2x", "4x"),
                         labsT=c("0.75x", "1x", "1.5x", "2x", "4x"),
                         titleX="hazard ratio",
                         titleT="relative risk of mandates",
                         pch=15, col="#E41A1C",
                         entryheight=0.19, topOffset=0.04, topIncr=0.0525,
                         plotwidth=1.25, filewidth=7,
                         sortBy=1, baseline=0, sims=10000, rev=FALSE,
                         logaxis=FALSE, robustness=FALSE, fancy=FALSE) {

    nTraces <- 1
    if (robustness) {
        if (is.list(xNum)) {
            nTraces <- length(xNum)
        }
        CF1 <- CF1ci90 <- vector("list", nTraces)
        for (i in 1:nTraces) {
            if (nTraces>1) {
                CF1[[i]] <- plotCovidCPHhelperRobust(res1, scenname, xNum[[i]], xpre1[[i]], xpost1[[i]],
                                                     baseline, sims=sims)
                CF1ci90[[i]] <- plotCovidCPHhelperRobust(res1, scenname, xNum[[i]], xpre1[[i]], xpost1[[i]],
                                                         baseline, ci=0.9, sims=sims)
            } else {
                CF1[[i]] <- plotCovidCPHhelperRobust(res1, scenname, xNum, xpre1, xpost1, baseline, sims=sims)
                CF1ci90[[i]] <- plotCovidCPHhelperRobust(res1, scenname, xNum, xpre1, xpost1,
                                                         baseline, ci=0.9, sims=sims)
            }
        }
    } else {
        CF1 <- plotCovidCPHhelper(res1, scenname, xNum, xpre1, xpost1, baseline, sims=sims)
        CF1ci90 <- plotCovidCPHhelper(res1, scenname, xNum, xpre1, xpost1, baseline, ci=0.9, sims=sims)
    }

    if (!is.null(res2)) {
        CF2 <- plotCovidCPHhelper(res2, scenname, xNum, xpre2, xpost2, baseline, sims=sims)
        CF2ci90 <- plotCovidCPHhelper(res2, scenname, xNum, xpre2, xpost2, baseline, ci=0.9, sims=sims)
    }

    ordering <- 1:length(scenname)

    if (!robustness) {
        if (sortBy==1) {
            ## Sort results by effect size for better plotting
            ordering <- order(CF1$CFs[,1])   
        }
        
        if (sortBy==2) {
            ## Sort results by effect size for better plotting
            ordering <- order(CF2$CFs[,1])    
        }

        if (rev) ordering <- rev(ordering)
    

        sortedCFnames <- CF1$CFnames[ordering]
        sortedCFs <- CF1$CFs[ordering,]
        sortedCFs90 <- CF1ci90$CFs[ordering,]

    }
    
    if (robustness) {
        shadowbox <- lighten(col, 0.65)
    } else {
        shadowbox <- NA
    }

    if (robustness) {

        trace1a <- traceSig1a <- traceSig1c <- vector("list", nTraces)
        
        for (j in 1:nTraces) {
            ## Main trace of simulation results
            curlabs <- CF1[[j]]$CFnames
            if (j>1) curlabs <- NULL
            curCol <- col[[j]]
            curPch <- pch[[j]]
            shadowbox <- lighten(curCol, 0.65)
            
            trace1a[[j]] <- ropeladder(x = CF1[[j]]$CFs[,1],
                                       lower = CF1[[j]]$CFs[,2],
                                       upper = CF1[[j]]$CFs[,3],
                                       labels = curlabs,
                                       entryheight=entryheight,  ## 0.19 .25
                                       ##subentryheight=0.7,  ## 1.1
                                       ##sublabels=sublabel1,
                                       ##sublabelsyoffset=0.025,  ## .0275
                                       pch=curPch,##c(15, 15, rep(16, length(sortedCFs)-2)),
                                       col=curCol,
                                       shadowbox=shadowbox,
                                       size=0.725,
                                       lex=1.75,
                                       lineend="square",
                                       group=1,
                                       plot=j
                                       )

            ## Mark statistically significant results
            sigMark1 <-  CF1[[j]]$CFs[,1]
            is.na(sigMark1) <- ((CF1[[j]]$CFs[,2]-1)*(CF1[[j]]$CFs[,3]-1)) > 0
            traceSig1a[[j]] <- ropeladder(x=sigMark1,
                                          col="white",
                                          pch=curPch,##c(15, 15, rep(16, length(sortedCFs)-1)),
                                          group=1,
                                          plot=j)
            
            ## Check for sig at 90% only; overplot a shaded circle
            sigMark1ci90 <- rep(NA, length(sigMark1))
            for (i in 1:length(sigMark1)) {
                if (!is.na(sigMark1[i])) {
                    if (((CF1ci90[[j]]$CFs[i,2]-1)*(CF1ci90[[j]]$CFs[i,3]-1)) > 0)
                        sigMark1ci90[i] <- CF1[[j]]$CFs[i,1]
                }
            }
            pcol <- rep(NA, length(curCol))
            for (i in 1:length(curCol)) pcol[i] <- lighten(curCol[i], pct=0.65)
            traceSig1c[[j]] <- ropeladder(x=sigMark1ci90,
                                          col=pcol,
                                          pch=curPch,##c(15, 15, rep(16, length(sortedCFs)-1)),
                                          layer=1,
                                          alpha=1,
                                          group=1,
                                          plot=j)
        }

    } else {
        ## Main trace of simulation results
        trace1a <- ropeladder(x = sortedCFs[,1],
                              lower = sortedCFs[,2],
                              upper = sortedCFs[,3],
                              labels = sortedCFnames,
                              entryheight=entryheight,  ## 0.19 .25
                              ##subentryheight=0.7,  ## 1.1
                              ##sublabels=sublabel1,
                              ##sublabelsyoffset=0.025,  ## .0275
                              pch=pch,##c(15, 15, rep(16, length(sortedCFs)-2)),
                              col=col,
                              shadowbox=shadowbox,
                              size=0.725,
                              lex=1.75,
                              lineend="square",
                              group=1,
                              plot=1
                              )
        
        ## Mark statistically significant results
        sigMark1 <- sortedCFs[,1]
        is.na(sigMark1) <- ((sortedCFs[,2]-1)*(sortedCFs[,3]-1)) > 0
        traceSig1a <- ropeladder(x=sigMark1,
                                 col="white",
                                 pch=pch,##c(15, 15, rep(16, length(sortedCFs)-1)),
                                 group=1,
                                 plot=1)
        
        
        ## Check for sig at 90% only; overplot a shaded circle
        sigMark1ci90 <- rep(NA, length(sigMark1))
        for (i in 1:length(sigMark1)) {
            if (!is.na(sigMark1[i])) {
                if (((sortedCFs90[i,2]-1)*(sortedCFs90[i,3]-1)) > 0)
                    sigMark1ci90[i] <- sortedCFs[i,1]
            }
        }
        pcol <- rep(NA, length(col))
        for (i in 1:length(col)) pcol[i] <- lighten(col[i], pct=0.65)
        traceSig1c <- ropeladder(x=sigMark1ci90,
                                 col=pcol,
                                 pch=pch,##c(15, 15, rep(16, length(sortedCFs)-1)),
                                 layer=1,
                                 alpha=1,
                                 group=1,
                                 plot=1)
        
    }

    collectedTopLabs <- NULL
    if (is.list(labsT)) {
        ctr <- 0
        collectedTopLabs <- vector("list", length(unlist(labsT)))
        for (i in 1:length(labsT)) {
            for (j in 1:length(labsT[[i]])) {
                ctr <- ctr+1
##                                                    y=1.035 + (j-1)*0.0475,

                
                collectedTopLabs[[ctr]] <- textTile(labels=rev(labsT[[i]])[j],
                                                    x=at[i],
                                                    y=1 + topOffset + (j-1)*topIncr,
                                                    clip="off",
                                                    cex=0.8,
                                                    plot=1:nTraces)
            }
        }
        labsT <- rep("", length(at))
    }
    
    ## Add reference line at 1.0
    vertmark <- linesTile(x=c(1,1), y=c(0,1), layer=12, plot=1:nTraces)

    if (fancy) {
        output <- list(outfile=file, width=filewidth, family="GillSans")
    } else {
        output <- list(outfile=file, width=filewidth) # 7
    }
    

         ## Make Figure 
        tile(trace1a, vertmark, traceSig1a,
         traceSig1c,
         collectedTopLabs,
         limits=limits,
         gridlines=list(type="xt"),
         topaxis=list(add=TRUE, log=logaxis, at=at,
                      labels=labsT),
         xaxis=list(log=logaxis, at=at,
                    labels=labsX),
         topaxistitle=list(labels=titleT, y=2.0),
         xaxistitle=list(labels=titleX, y=0.5),
         maintitle=list(labels=""),
         width=list(plot=plotwidth),
         height=list(plottitle=3,xaxistitle=3.0,topaxistitle=3.0),
         output=output
         )

  

    if (fancy) embedGS(paste(file,".pdf",sep=""))

    if (!robustness) {
        list(scen=sortedCFnames, hr95=sortedCFs, hr90=sortedCFs90)
    }
}

plotCovidCPHhelper <- function(res, scenname, xNum, xpre, xpost,
                               baseline=0, ci=0.95, sims=10000) {
    
    ## Simulate coefficients to account for estimation uncertainty
    simbetas <- exp(mvrnorm(sims, coef(res), vcov(res)))
    
    ## Collect counterfactual changes for each covariate:
    ## typically IQR range for continuous covariates
    ## and 0 to 1 change for binary covariates
    CFnames <- CFs <- NULL

    for (i in 1:length(scenname)) {
        CFnames <- c(CFnames, scenname[i])
        CFs <- rbind(CFs, simCoxPH(simbetas[,xNum[i]], xpre[i], xpost[i], 0, ci=ci))
    }

    list(CFnames=CFnames, CFs=CFs)
}



plotCovidCPHhelperRobust <- function(res, scenname, xNum, xpre, xpost,
                                     baseline=0, ci=0.95, sims=10000) {

    ## Collect counterfactual changes for each covariate:
    ## typically IQR range for continuous covariates
    ## and 0 to 1 change for binary covariates
    CFnames <- CFs <- NULL
    
    nTests <- length(res)
    for (i in 1:nTests) {

        ## check for alternative scenarios
        if (length(xpre)>1) {
            curPre <- xpre[i]
            curPost <- xpost[i]
        } else {
            curPre <- xpre
            curPost <- xpost
        }

        ## check for varying location of coefficient of interest
        if (length(xNum)>1) {
            curxNum <- xNum[i]
        } else {
            curxNum <- xNum
        }

        ## Simulate coefficients to account for estimation uncertainty
        simbetas <- exp(mvrnorm(sims, coef(res[[i]]), vcov(res[[i]])))
    
        CFnames <- c(CFnames, scenname[i])
        if (is.na(curxNum)) {
            CFs <- rbind(CFs, rep(NA, 3))
        } else {
            CFs <- rbind(CFs, simCoxPH(simbetas[,curxNum], curPre, curPost, 0, ci=ci))
        }

    }
    
    list(CFnames=CFnames, CFs=CFs)
}



## Plot average marginal effects from Cox PH models using coxed and tile
coxphAMEplot <- function(file, model, data, id, nameScen, cf1, cf2, sims=200, method="npsf", ciMethod="bca",
                         rev=TRUE, pch=pch, col="#E41A1C", log=TRUE,
                         limits=c(5/24, 3.25),
                         at=c(1/4, 0.5, 1, 2, 3),
                         labs=c("6 hours", "12 hours", "1 day", "2 days", "3 days"),
                         title="Additional Expected Delay in Mandating Social Distancing",
                         maxit=15,
                         fancy=FALSE) {

    
    mdata <- extractdata(model, data, extra=NULL, na.rm = TRUE)
    nScen <- length(cf1)
    me95ALL  <- me90ALL <- NULL
    for (i in 1:nScen) {

        eval(str2lang(
            paste0("me95 <- coxed2(model, id=id, method=method, bootstrap = TRUE, B = sims, level = .95, maxit=maxit, confidence=ciMethod, newdata = dplyr::mutate(mdata,",
                   cf1[i],
                   "), newdata2 = dplyr::mutate(mdata,",
                   cf2[i],
                   "))")))

        
        eval(str2lang(
            paste0("me90 <- coxed2(model, id=id, method=method, bootstrap = TRUE, B = sims, level = .90, maxit=maxit, confidence=ciMethod, newdata = dplyr::mutate(mdata,",
                   cf1[i],
                   "), newdata2 = dplyr::mutate(mdata,",
                   cf2[i],
                   "))")))
        
        me95ALL <- rbind(me95ALL,
                         c(me95$mean.diff$mean, me95$mean.diff$lb, me95$mean.diff$ub))
        me90ALL <- rbind(me95ALL,
                       c(me90$mean.diff$mean, me90$mean.diff$lb, me90$mean.diff$ub))
    }

    ##list(me95=me95ALL, me90=me90ALL)

    ## Convert from duration to delay
    me95ALL <- -me95ALL
    me90ALL <- -me90ALL
    
    ## Sort results by effect size for better plotting
    ordering <- order(me95ALL[,1])   

    if (rev) ordering <- rev(ordering)

    sortedCFnames <- nameScen[ordering]
    sortedCFs <- me95ALL[ordering,]
    sortedCFs90 <- me90ALL[ordering,]

    
    ## Main trace of simulation results
    trace1a <- ropeladder(x = sortedCFs[,1],
                          lower = sortedCFs[,2],
                          upper = sortedCFs[,3],
                          labels = sortedCFnames,
                          entryheight=0.19,  ## .25
                          ##subentryheight=0.7,  ## 1.1
                          #sublabels=sublabel1,
                          #sublabelsyoffset=0.025,  ## .0275
                          pch=pch,##c(15, 15, rep(16, length(sortedCFs)-2)),
                          col=col,
                          size=0.725,
                          lex=1.75,
                          lineend="square",
                          group=1,
                          plot=1
                          )           
    
    
    ## Mark statistically significant results
    sigMark1 <- sortedCFs[,1]
    is.na(sigMark1) <- ((sortedCFs[,2])*(sortedCFs[,3])) > 0
    traceSig1a <- ropeladder(x=sigMark1,
                             col="white",
                             pch=pch,##c(15, 15, rep(16, length(sortedCFs)-1)),
                             group=1,
                             plot=1)

   
    ## Check for sig at 90% only; overplot a shaded circle
    sigMark1ci90 <- rep(NA, length(sigMark1))
    for (i in 1:length(sigMark1)) {
        if (!is.na(sigMark1[i])) {
            if (((sortedCFs90[i,2])*(sortedCFs90[i,3])) > 0)
                sigMark1ci90[i] <- sortedCFs[i,1]
        }
    }
    pcol <- rep(NA, length(col))
    for (i in 1:length(col)) pcol[i] <- lighten(col[i], pct=0.65)
    traceSig1c <- ropeladder(x=sigMark1ci90,
                             col=pcol,
                             pch=pch,##c(15, 15, rep(16, length(sortedCFs)-1)),
                             layer=1,
                             alpha=1,
                             group=1,
                             plot=1)
    
    collectedTopLabs <- NULL
    if (is.list(labs)) {
        ctr <- 0
        collectedTopLabs <- vector("list", length(unlist(labs)))
        for (i in 1:length(labs)) {
            for (j in 1:length(labs[[i]])) {
                ctr <- ctr+1
                collectedTopLabs[[ctr]] <- textTile(labels=rev(labs[[i]])[j],
                                                    x=at[i],
                                                    y=1.035 + (j-1)*0.0475,
                                                    clip="off",
                                                    cex=0.8,
                                                    plot=1)
            }
        }
        labs <- rep("", length(at))
    }
    

    ## Add reference line at .0
    vertmark <- linesTile(x=c(0,0), y=c(0,1), layer=12, plot=1)


    if (fancy) {
        output <- list(outfile=file, width=7, family="GillSans")
    } else {
        output <- list(outfile=file, width=7)
    }
    
    ## Make Figure 
    tile(trace1a, #vertmark,
         traceSig1a,
         traceSig1c,
         collectedTopLabs,
         limits=limits,
         gridlines=list(type="xt"),
         topaxis=list(add=TRUE, log=log, at=at,
                      labels=labs),
         xaxis=list(log=log, at=at,
                   labels=rep("", length(at))),
         topaxistitle=list(labels=title, y=1.75),
         ##xaxistitle=list(labels=titleX, y=0.5),
         ##maintitle=list(labels=""),
         width=list(plot=1.25),
         height=list(plottitle=3,xaxistitle=3.0,topaxistitle=3.0),
         output=output
         )
    if (fancy) embedGS(paste(file,".pdf",sep=""))

    list(me95=me95ALL, me90=me90ALL)
}



#######################################################################
## When originally written, this paper used version 0.3.2 of the coxed package,
## which began to produce errors when coxed updated to 0.3.3.
## 
## This (current) version of coxed appears to have a bug when the 
## bootstrap option with newdata.  In particular coxed() crashes on this line:
##
##          boot.cph <- rms::cph(S ~ X, x = TRUE, y = TRUE)
##
## It appears this line should be instead:
##
##          boot.cph <- rms::cph(S ~ Xmat, x = TRUE, y = TRUE)
##
## The authors provide the replacement function coxed2 below to correct this.
## The place where the line has been repaired is noted; otherwise the function
## is identical to coxed in the coxed package version 0.3.3.
##
## The coxed package maintainers have been notified.

coxed2 <- function (cox.model, newdata = NULL, newdata2 = NULL, bootstrap = FALSE, 
    method = "npsf", k = -1, B = 200, confidence = "studentized", 
    level = 0.95, id = NULL, ...) 
{
    tvc <- (ncol(cox.model$y) == 3)
    if (is.null(newdata) & !is.null(newdata2)) 
        stop("newdata2 must only be specified if newdata is also specified")
    if (!method %in% c("gam", "npsf")) 
        stop("method must be one of 'gam' and 'npsf'")
    if (!is.null(newdata2) & !all(dim(newdata) == dim(newdata2))) 
        stop("newdata and newdata2 must have the same dimensions")
    if (!confidence %in% c("studentized", "empirical", "bca")) 
        stop("confidence must be one of 'studentized', 'empirical', or 'bca'")
    if (tvc & is.null(id)) 
        stop("id must be filled in with the case ID in the data used to estimate the Cox PH model if you have time-varying covariates.")
    if (method == "gam") {
        if (!tvc) 
            dur1 <- coxed.gam(cox.model, newdata = newdata, k = k)
        if (tvc) 
            dur1 <- coxed.gam.tvc(cox.model, newdata = newdata, 
                k = k)
        dur1.pe <- dur1$exp.dur
        gam.model <- dur1$gam.model
        gam.data <- dur1$gam.data
        m1 <- mean(gam.data$y, na.rm = TRUE)
        m2 <- mean(gam.data$y[gam.data$failed == 1], na.rm = TRUE)
        if (m1/m2 > 1.25) 
            warning("Inclusion of censored cases increases mean duration by more than 25%. Consider using method='npsf' instead")
    }
    if (method == "npsf") {
        if (!tvc) 
            dur1 <- coxed.npsf(cox.model, newdata = newdata)
        if (tvc) 
            dur1 <- coxed.npsf.tvc(cox.model, newdata = newdata)
        dur1.pe <- dur1$exp.dur
        baseline.functions <- dur1$baseline.functions
    }
    exp.dur <- data.frame(exp.dur = dur1.pe)
    if (!is.null(newdata2)) {
        if (method == "gam" & !tvc) 
            dur2.pe <- coxed.gam(cox.model, newdata = newdata2, 
                k = k)$exp.dur
        if (method == "gam" & tvc) 
            dur2.pe <- coxed.gam.tvc(cox.model, newdata = newdata2, 
                k = k)$exp.dur
        if (method == "npsf" & !tvc) 
            dur2.pe <- coxed.npsf(cox.model, newdata = newdata2)$exp.dur
        if (method == "npsf" & tvc) 
            dur2.pe <- coxed.npsf.tvc(cox.model, newdata = newdata2)$exp.dur
        diff <- dur2.pe - dur1.pe
        exp.dur <- data.frame(exp.dur1 = dur1.pe, exp.dur2 = dur2.pe, 
            difference = diff)
    }
    if (!bootstrap) {
        res <- list(exp.dur = exp.dur, mean = colMeans(exp.dur, 
            na.rm = TRUE), median = apply(exp.dur, 2, FUN = function(x) {
            median(x, na.rm = TRUE)
        }))
    }
    if (bootstrap) {
        if (!tvc) {
            S <- Surv(as.numeric(cox.model$y[, 1]), as.numeric(cox.model$y[, 
                2]), type = "right")
            Xmat <- model.matrix(cox.model)
            boot.model <- bootcov2(rms::cph(S ~ Xmat, x = TRUE, 
                y = TRUE), B = B, ...)
        }
        if (tvc) {
            id <- id[as.numeric(rownames(model.frame(cox.model)))]
            S <- Surv(as.numeric(cox.model$y[, 1]), as.numeric(cox.model$y[, 
                2]), as.numeric(cox.model$y[, 3]), type = "counting")
            Xmat <- model.matrix(cox.model)

######################################################
## This line repaired:  originally had X not Xmat
            boot.cph <- rms::cph(S ~ Xmat, x = TRUE, y = TRUE)
######################################################
            class(boot.cph) <- "tvc"
            boot.model <- bootcov2(boot.cph, B = B, cluster = id, 
                                   ...)

        }
        bs.coef <- boot.model$boot.Coef
        bs.obs <- na.omit(boot.model$b.ind)
        if (is.null(newdata2)) {
            exp.dur.mat <- matrix(numeric(), nrow(exp.dur), nrow(bs.coef))
            for (i in 1:nrow(bs.coef)) {
                if (method == "gam" & !tvc) 
                  dur1.pe <- coxed.gam(cox.model, newdata = newdata, 
                    warn = FALSE, k = k, coef = bs.coef[i, ], 
                    b.ind = bs.obs[, i])$exp.dur
                if (method == "gam" & tvc) 
                  dur1.pe <- coxed.gam.tvc(cox.model, newdata = newdata, 
                    warn = FALSE, k = k, coef = bs.coef[i, ], 
                    b.ind = bs.obs[, i])$exp.dur
                if (method == "npsf" & !tvc) 
                  dur1.pe <- coxed.npsf(cox.model, newdata = newdata, 
                    coef = bs.coef[i, ], b.ind = bs.obs[, i])$exp.dur
                if (method == "npsf" & tvc) 
                  dur1.pe <- coxed.npsf.tvc(cox.model, newdata = newdata, 
                    coef = bs.coef[i, ], b.ind = bs.obs[, i])$exp.dur
                exp.dur.mat[, i] <- dur1.pe
            }
            exp.dur.se <- apply(exp.dur.mat, 1, FUN = function(x) {
                sd(x, na.rm = TRUE)
            })
            mean.vec <- colMeans(exp.dur.mat, na.rm = TRUE)
            mean.se <- sd(mean.vec, na.rm = TRUE)
            median.vec <- apply(exp.dur.mat, 2, FUN = function(x) {
                median(x, na.rm = TRUE)
            })
            median.se <- sd(median.vec, na.rm = TRUE)
            if (confidence == "empirical") {
                exp.dur.lb <- apply(exp.dur.mat, 1, FUN = function(x) {
                  quantile(x, (1 - level)/2, names = FALSE, na.rm = TRUE)
                })
                exp.dur.ub <- apply(exp.dur.mat, 1, FUN = function(x) {
                  quantile(x, level + (1 - level)/2, names = FALSE, 
                    na.rm = TRUE)
                })
                mean.lb <- quantile(mean.vec, (1 - level)/2, 
                  names = FALSE, na.rm = TRUE)
                mean.ub <- quantile(mean.vec, level + (1 - level)/2, 
                  names = FALSE, na.rm = TRUE)
                median.lb <- quantile(median.vec, (1 - level)/2, 
                  names = FALSE, na.rm = TRUE)
                median.ub <- quantile(median.vec, level + (1 - 
                  level)/2, names = FALSE, na.rm = TRUE)
            }
            else if (confidence == "studentized") {
                exp.dur.lb <- exp.dur$exp.dur + qnorm((1 - level)/2) * 
                  exp.dur.se
                exp.dur.ub <- exp.dur$exp.dur + qnorm(level + 
                  (1 - level)/2) * exp.dur.se
                mean.lb <- mean(exp.dur$exp.dur, na.rm = TRUE) + 
                  qnorm((1 - level)/2) * mean.se
                mean.ub <- mean(exp.dur$exp.dur, na.rm = TRUE) + 
                  qnorm(level + (1 - level)/2) * mean.se
                median.lb <- median(exp.dur$exp.dur, na.rm = TRUE) + 
                  qnorm((1 - level)/2) * median.se
                median.ub <- median(exp.dur$exp.dur, na.rm = TRUE) + 
                  qnorm(level + (1 - level)/2) * median.se
            }
            else if (confidence == "bca") {
                exp.dur.bca <- apply(exp.dur.mat, 1, FUN = function(x) {
                  bca(x, conf.level = level)
                })
                exp.dur.lb <- exp.dur.bca[1, ]
                exp.dur.ub <- exp.dur.bca[2, ]
                mean.bca <- bca(mean.vec, conf.level = level)
                mean.lb <- mean.bca[1]
                mean.ub <- mean.bca[2]
                median.bca <- bca(median.vec, conf.level = level)
                median.lb <- median.bca[1]
                median.ub <- median.bca[2]
            }
            res <- list(exp.dur = data.frame(exp.dur = exp.dur$exp.dur, 
                bootstrap.se = exp.dur.se, lb = exp.dur.lb, ub = exp.dur.ub), 
                mean = data.frame(mean = mean(exp.dur$exp.dur, 
                  na.rm = TRUE), bootstrap.se = mean.se, lb = mean.lb, 
                  ub = mean.ub), median = data.frame(median = median(exp.dur$exp.dur, 
                  na.rm = TRUE), bootstrap.se = median.se, lb = median.lb, 
                  ub = median.ub))
        }
        else {
            exp.dur1.mat <- matrix(numeric(), nrow(exp.dur), 
                nrow(bs.coef))
            exp.dur2.mat <- matrix(numeric(), nrow(exp.dur), 
                nrow(bs.coef))
            diff.mat <- matrix(numeric(), nrow(exp.dur), nrow(bs.coef))
            for (i in 1:nrow(bs.coef)) {
                if (method == "gam" & !tvc) {
                  warn <- (i == 1)
                  dur1.pe <- coxed.gam(cox.model, newdata = newdata, 
                    k = k, warn = warn, coef = bs.coef[i, ], 
                    b.ind = bs.obs[, i])$exp.dur
                  dur2.pe <- coxed.gam(cox.model, newdata = newdata2, 
                    k = k, warn = warn, coef = bs.coef[i, ], 
                    b.ind = bs.obs[, i])$exp.dur
                }
                if (method == "gam" & tvc) {
                  warn <- (i == 1)
                  dur1.pe <- coxed.gam.tvc(cox.model, newdata = newdata, 
                    k = k, warn = warn, coef = bs.coef[i, ], 
                    b.ind = bs.obs[, i])$exp.dur
                  dur2.pe <- coxed.gam.tvc(cox.model, newdata = newdata2, 
                    k = k, warn = warn, coef = bs.coef[i, ], 
                    b.ind = bs.obs[, i])$exp.dur
                }
                if (method == "npsf" & !tvc) {
                  dur1.pe <- coxed.npsf(cox.model, newdata = newdata, 
                    coef = bs.coef[i, ], b.ind = bs.obs[, i])$exp.dur
                  dur2.pe <- coxed.npsf(cox.model, newdata = newdata2, 
                    coef = bs.coef[i, ], b.ind = bs.obs[, i])$exp.dur
                }
                if (method == "npsf" & tvc) {
                  dur1.pe <- coxed.npsf.tvc(cox.model, newdata = newdata, 
                    coef = bs.coef[i, ], b.ind = bs.obs[, i])$exp.dur
                  dur2.pe <- coxed.npsf.tvc(cox.model, newdata = newdata2, 
                    coef = bs.coef[i, ], b.ind = bs.obs[, i])$exp.dur
                }
                exp.dur1.mat[, i] <- dur1.pe
                exp.dur2.mat[, i] <- dur2.pe
                diff.mat[, i] <- dur2.pe - dur1.pe
            }
            exp.dur1.se <- apply(exp.dur1.mat, 1, FUN = function(x) {
                sd(x, na.rm = TRUE)
            })
            mean1.vec <- colMeans(exp.dur1.mat, na.rm = TRUE)
            mean1.se <- sd(mean1.vec, na.rm = TRUE)
            median1.vec <- apply(exp.dur1.mat, 2, FUN = function(x) {
                median(x, na.rm = TRUE)
            })
            median1.se <- sd(median1.vec, na.rm = TRUE)
            exp.dur2.se <- apply(exp.dur2.mat, 1, FUN = function(x) {
                sd(x, na.rm = TRUE)
            })
            mean2.vec <- colMeans(exp.dur2.mat, na.rm = TRUE)
            mean2.se <- sd(mean2.vec, na.rm = TRUE)
            median2.vec <- apply(exp.dur2.mat, 2, FUN = function(x) {
                median(x, na.rm = TRUE)
            })
            median2.se <- sd(median2.vec, na.rm = TRUE)
            diff.se <- apply(diff.mat, 1, FUN = function(x) {
                sd(x, na.rm = TRUE)
            })
            mean.diff.vec <- colMeans(diff.mat, na.rm = TRUE)
            mean.diff.se <- sd(mean.diff.vec, na.rm = TRUE)
            median.diff.vec <- apply(diff.mat, 2, FUN = function(x) {
                median(x, na.rm = TRUE)
            })
            median.diff.se <- sd(median.diff.vec, na.rm = TRUE)
            if (confidence == "empirical") {
                exp.dur1.lb <- apply(exp.dur1.mat, 1, FUN = function(x) {
                  quantile(x, (1 - level)/2, names = FALSE, na.rm = TRUE)
                })
                exp.dur1.ub <- apply(exp.dur1.mat, 1, FUN = function(x) {
                  quantile(x, level + (1 - level)/2, names = FALSE, 
                    na.rm = TRUE)
                })
                mean1.lb <- quantile(mean1.vec, (1 - level)/2, 
                  names = FALSE, na.rm = TRUE)
                mean1.ub <- quantile(mean1.vec, level + (1 - 
                  level)/2, names = FALSE, na.rm = TRUE)
                median1.lb <- quantile(median1.vec, (1 - level)/2, 
                  names = FALSE, na.rm = TRUE)
                median1.ub <- quantile(median1.vec, level + (1 - 
                  level)/2, names = FALSE, na.rm = TRUE)
                exp.dur2.lb <- apply(exp.dur2.mat, 1, FUN = function(x) {
                  quantile(x, (1 - level)/2, names = FALSE, na.rm = TRUE)
                })
                exp.dur2.ub <- apply(exp.dur2.mat, 1, FUN = function(x) {
                  quantile(x, level + (1 - level)/2, names = FALSE, 
                    na.rm = TRUE)
                })
                mean2.lb <- quantile(mean2.vec, (1 - level)/2, 
                  names = FALSE, na.rm = TRUE)
                mean2.ub <- quantile(mean2.vec, level + (1 - 
                  level)/2, names = FALSE, na.rm = TRUE)
                median2.lb <- quantile(median2.vec, (1 - level)/2, 
                  names = FALSE, na.rm = TRUE)
                median2.ub <- quantile(median2.vec, level + (1 - 
                  level)/2, names = FALSE, na.rm = TRUE)
                diff.lb <- apply(diff.mat, 1, FUN = function(x) {
                  quantile(x, (1 - level)/2, names = FALSE, na.rm = TRUE)
                })
                diff.ub <- apply(diff.mat, 1, FUN = function(x) {
                  quantile(x, level + (1 - level)/2, names = FALSE, 
                    na.rm = TRUE)
                })
                mean.diff.lb <- quantile(mean.diff.vec, (1 - 
                  level)/2, names = FALSE, na.rm = TRUE)
                mean.diff.ub <- quantile(mean.diff.vec, level + 
                  (1 - level)/2, names = FALSE, na.rm = TRUE)
                median.diff.lb <- quantile(median.diff.vec, (1 - 
                  level)/2, names = FALSE, na.rm = TRUE)
                median.diff.ub <- quantile(median.diff.vec, level + 
                  (1 - level)/2, names = FALSE, na.rm = TRUE)
            }
            else if (confidence == "studentized") {
                exp.dur1.lb <- exp.dur[, 1] + qnorm((1 - level)/2) * 
                  exp.dur1.se
                exp.dur1.ub <- exp.dur[, 1] + qnorm(level + (1 - 
                  level)/2) * exp.dur1.se
                mean1.lb <- mean(exp.dur[, 1], na.rm = TRUE) + 
                  qnorm((1 - level)/2) * mean1.se
                mean1.ub <- mean(exp.dur[, 1], na.rm = TRUE) + 
                  qnorm(level + (1 - level)/2) * mean1.se
                median1.lb <- median(exp.dur[, 1], na.rm = TRUE) + 
                  qnorm((1 - level)/2) * median1.se
                median1.ub <- median(exp.dur[, 1], na.rm = TRUE) + 
                  qnorm(level + (1 - level)/2) * median1.se
                exp.dur2.lb <- exp.dur[, 2] + qnorm((1 - level)/2) * 
                  exp.dur2.se
                exp.dur2.ub <- exp.dur[, 2] + qnorm(level + (1 - 
                  level)/2) * exp.dur2.se
                mean2.lb <- mean(exp.dur[, 2], na.rm = TRUE) + 
                  qnorm((1 - level)/2) * mean2.se
                mean2.ub <- mean(exp.dur[, 2], na.rm = TRUE) + 
                  qnorm(level + (1 - level)/2) * mean2.se
                median2.lb <- median(exp.dur[, 2], na.rm = TRUE) + 
                  qnorm((1 - level)/2) * median2.se
                median2.ub <- median(exp.dur[, 2], na.rm = TRUE) + 
                  qnorm(level + (1 - level)/2) * median2.se
                diff.lb <- exp.dur[, 3] + qnorm((1 - level)/2) * 
                  diff.se
                diff.ub <- exp.dur[, 3] + qnorm(level + (1 - 
                  level)/2) * diff.se
                mean.diff.lb <- mean(exp.dur[, 3], na.rm = TRUE) + 
                  qnorm((1 - level)/2) * mean.diff.se
                mean.diff.ub <- mean(exp.dur[, 3], na.rm = TRUE) + 
                  qnorm(level + (1 - level)/2) * mean.diff.se
                median.diff.lb <- median(exp.dur[, 3], na.rm = TRUE) + 
                  qnorm((1 - level)/2) * median.diff.se
                median.diff.ub <- median(exp.dur[, 3], na.rm = TRUE) + 
                  qnorm(level + (1 - level)/2) * median.diff.se
            }
            else if (confidence == "bca") {
                exp.dur1.bca <- apply(exp.dur1.mat, 1, FUN = function(x) {
                  bca(x, conf.level = level)
                })
                exp.dur1.lb <- exp.dur1.bca[1, ]
                exp.dur1.ub <- exp.dur1.bca[2, ]
                mean1.bca <- bca(mean1.vec, conf.level = level)
                mean1.lb <- mean1.bca[1]
                mean1.ub <- mean1.bca[2]
                median1.bca <- bca(median1.vec, conf.level = level)
                median1.lb <- median1.bca[1]
                median1.ub <- median1.bca[2]
                exp.dur2.bca <- apply(exp.dur2.mat, 1, FUN = function(x) {
                  bca(x, conf.level = level)
                })
                exp.dur2.lb <- exp.dur2.bca[1, ]
                exp.dur2.ub <- exp.dur2.bca[2, ]
                mean2.bca <- bca(mean2.vec, conf.level = level)
                mean2.lb <- mean2.bca[1]
                mean2.ub <- mean2.bca[2]
                median2.bca <- bca(median2.vec, conf.level = level)
                median2.lb <- median2.bca[1]
                median2.ub <- median2.bca[2]
                diff.bca <- apply(diff.mat, 1, FUN = function(x) {
                  bca(x, conf.level = level)
                })
                diff.lb <- diff.bca[1, ]
                diff.ub <- diff.bca[2, ]
                mean.diff.bca <- bca(mean.diff.vec, conf.level = level)
                mean.diff.lb <- mean.diff.bca[1]
                mean.diff.ub <- mean.diff.bca[2]
                median.diff.bca <- bca(median.diff.vec, conf.level = level)
                median.diff.lb <- median.diff.bca[1]
                median.diff.ub <- median.diff.bca[2]
            }
            res <- list(exp.dur1 = data.frame(exp.dur = exp.dur[, 
                1], bootstrap.se = exp.dur1.se, lb = exp.dur1.lb, 
                ub = exp.dur1.ub), mean1 = data.frame(mean = mean(exp.dur[, 
                1], na.rm = TRUE), bootstrap.se = mean1.se, lb = mean1.lb, 
                ub = mean1.ub), median1 = data.frame(median = median(exp.dur[, 
                1], na.rm = TRUE), bootstrap.se = median1.se, 
                lb = median1.lb, ub = median1.ub), exp.dur2 = data.frame(exp.dur = exp.dur[, 
                2], bootstrap.se = exp.dur2.se, lb = exp.dur2.lb, 
                ub = exp.dur2.ub), mean2 = data.frame(mean = mean(exp.dur[, 
                2], na.rm = TRUE), bootstrap.se = mean2.se, lb = mean2.lb, 
                ub = mean2.ub), median2 = data.frame(median = median(exp.dur[, 
                2], na.rm = TRUE), bootstrap.se = median2.se, 
                lb = median2.lb, ub = median2.ub), diff = data.frame(exp.dur = exp.dur[, 
                3], bootstrap.se = diff.se, lb = diff.lb, ub = diff.ub), 
                mean.diff = data.frame(mean = mean(exp.dur[, 
                  3], na.rm = TRUE), bootstrap.se = mean.diff.se, 
                  lb = mean.diff.lb, ub = mean.diff.ub), median.diff = data.frame(median = median(exp.dur[, 
                  3], na.rm = TRUE), bootstrap.se = median.diff.se, 
                  lb = median.diff.lb, ub = median.diff.ub))
        }
    }
    if (method == "gam") {
        res$gam.model <- gam.model
        res$gam.data <- gam.data
    }
    if (method == "npsf") {
        res$baseline.functions <- baseline.functions
    }
    if (!is.null(newdata2)) {
        class(res) <- "coxedMargin"
    }
    else {
        class(res) <- "coxedExpdur"
    }
    return(res)
}





